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SUMMARY 


The numerical model for a rocket thermal analysis code (RTE) is discussed. RTE is a 
comprehensive thermal analysis code for thermal analysis of regeneratively cooled rocket en- 
gines. The input to the code consists of the composition of fuel/oxidant mixture and flow rates, 
chamber pressure, coolant temperature and pressure, dimensions of the engine, materials and 
the number of nodes in different parts of the engine. The code allows for temperature variation 
in axial, radial and circumferential directions. By implementing an iterative scheme, it pro- 
vides nodal temperature distribution, rates of heat transfer, hot gas and coolant thermal and 
transport properties. The fuel/oxidant mixture ratio can be varied along the thrust chamber. 
This feature allows the user to incorporate a non-equilibrium model or an energy release model 
for the hot-gas-side. The user has the option of bypassing the hot-gas-side calculations and 
directly inputting the gas-side fluxes. This feature is used to link RTE to a boundary layer 
module for the hot-gas-side heat flux calculations. 


INTRODUCTION 

Thermal analysis is an essential and integral part in the design of rocket engines. The 
need for thermal analysis is especially important in the reusable engines where an effective 
and efficient cooling system becomes a crucial factor in extending the engine life. In the 
new high pressure engines, such as chemical transfer vehicle engines, hot-gas temperature is 
very high (can reach 7000R at the throat). It is therefore essential to be able to estimate 
the wall temperature and ensure that the material can withstand such high temperature. 
Furthermore, an accurate thermal model enables an engine designer to modify the cooling 
channel configuration for the maximum cooling at high temperature areas. 

The thermal phenomena in rocket engines involve interactions among a number of processes, 
including, combustion in the thrust chamber, expansion of hot-gases through the nozzle, heat 
transfer from hot-gases to the nozzle wall via convection and radiation, conduction in the 
wall, and convection to the cooling channel. Further complexities of the thermal analysis 
in rocket engines are due to three-dimensional geometry, coolant and hot gas heat transfer 
coefficient dependence on the pressure and wall temperature, unknown coolant pressure drop 
and properties, axial conduction of heat within the wall, and radiative heat transfer between 
gases and surfaces of the engine. A comprehensive thermal model must account for all of these 
items. 

RTE [1] is a comprehensive rocket thermal analysis code that uses a number of existing 
codes and allows interaction among them via some iterative procedures. The code is based 
on the geometry of a typical regeneratively-cooled engine similar to that shown in Figure 
h It uses CET (Chemical Equilibrium with Transport Properties) [2] and GASP [3] for the 
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evaluation of hot-gas and coolant properties. The inputs to this code consist of the composition 
of fuel/oxidant mixtures and flow rates, chamber pressure, coolant entrance temperature and 
pressure, dimensions of the engine and materials in different parts of the engine, as well as 
the mesh generation data. This program allows temperature variations in axial, radial and 
circumferential directions, and by implementing an iterative scheme it provides temperature 
distributions, rates of heat transfer, and hot-gas and coolant thermal and transport properties. 
The fuel/oxidant mixture ratio can be varied along the thrust chamber. This feature allows the 
user to incorporate a nonequilibrium model or an energy release model for the hot-gas-side. The 
mixture ratio along the thrust chamber is calculated using ROCCID [4] (ROCket Combustor 
Interactive Design and Analysis Computer Program). ROCCID has been modified to take RTE 
input and make the mixture ratio variable along the thrust chamber. The user has the option 
of bypassing the hot-gas-side calculations and directly inputting gas side fluxes. This feature 
is used to link RTE to a boundary layer program for the hot-gas-side heat flux calculation. 
The procedure for linking RTE to a hot-gas side program. TDK [5] (Two-Dimensional Kinetics 
Nozzle Performance Computer Program) is described here. 




Figure 1. A Rocket Thrust Chamber and Nozzle 


NUMERICAL MODEL 

The numerical model of the RTE is based on the geometry of a typical regeneratively- 
cooled thrust chamber (shown in Figure 1). The wall can consist of three layers: a coating, the 
channel, and the closeout. These three layers can be different materials or the same material. 
The number of cooling channels in the wall are also specified by the user. For the numerical 
procedure, the rocket thrust chamber and nozzle are subdivided into a number of stations 
along the longitudinal direction, as shown in Figure 2. The thermodynamic and transport 
properties of the combustion gases are evaluated using the chemical equilibrium composition 
computer program developed by Gordon and McBride [2. 6] (CET, Chemical Equilibrium 
with Transport properties). The GASP (GAS Properties) [3] or WASP (Water And Steam 
Properties) [7] programs are implemented to obtain coolant thermodynamic and transport 
properties. Since the heat transfer coefficients of the hot-gas and coolant sides are related to 
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surface temperatures, an iterative procedure is used to evaluate heat transfer coefficients and 
adiabatic wall temperatures. 


STATION NUMBERS 



Figure 2. A Rocket Thrust Chamber Subdivided into a Number of Stations 

The temperature distribution within the wall is determined via a three-dimensional finite 
difference scheme. In this method, finite difference grids are superimposed throughout the wall 
at different stations. The temperature of each node is then written in terms of temperatures of 
neighboring nodes (the four closest nodes at the same station and two nodes at the neighboring 
stations). The program marches axially from one station to another. At each station the Gauss- 
Siedel iterative method is used to obtain convergence for the temperature distribution along 
the radial and circumferential directions. When the axial march is completed, comparison 
is made between the results of the present march and that of the previous one to see if the 
convergence criteria in the axial direction has been met. If it is not met. the code starts again 
at the first station and makes another axial march. The process continues until convergence 
is achieved. A detailed description of this numerical model is outlined below. 

First, the static pressures, temperatures, enthalpies and Mach numbers for the combustion 
gases are evaluated using the ROCKET subroutine from [2]. It should be noted that these 
properties are independent of wall temperature and only depend on the cross-sectional area 
of the nozzle, the propellant used and chamber pressure. Indeed, the heat transfer from hot 
gases to the chamber and nozzle wall will cause very little change in the gas temperature (the 
thermodynamic process dominates the transport process). 

On the coolant side, the stagnation enthalpy and density at the entrance to the cooling 
channel are evaluated as functions of the coolant stagnation pressure and temperature (ico — 
ico(Pco, T C o) and pco = PCo(Pcq> Tco)) using the GASP or WASP programs. 

The model now begins its axial marches (passes) starting from the first station. At the 
first axial march an initial guess for the wall temperature distribution is made. For the next 
march, however, the results of temperature distribution for the previous march can be used as 
an initial guess. The hot gas and coolant adiabatic wall temperatures and wall properties can 
be evaluated at a given station based on the assumed wall temperature distribution using the 
properties computer codes [2, 6, 3, 7] for the combustion gases and the coolant. The reference 
enthalpy of the gas side. icx n is g iven b Y [8] 
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iGXn = °- 5 (tGW'„ + *GS n ) + 0 . 180 (ico n ~ iGSn ) (!) 

where iGW n is a function of gas static pressure Pcs„ and gas-side wall temperature Tew n and 
is evaluated using the program given in [2]. The gas-side adiabatic wall enthalpy, iGAW n is 
calculated using the following equation [8, 9] 


iGAW„ = iGSn + ( Pr GX n ) 1/3 (iG0„ ~ *GS„ ) 
where the gas reference Prandtl number Prcx n i s 


( 2 ) 


P^GXn 


C'PGXnV'GX 

k GXn 


(3) 


C pcx fJ-GXn a nd kGXn are functions of Pcs n an d iGXn- Once the gas-side adiabatic wall 
temperature is determined, the wall adiabatic temperature is calculated via 


TGAWn = f( P GS n ,iGAW n ) (4) 

and using the combustion codes 12. 6). The hot-gas side heat transfer coefficient. hc n is given 

by [8] 


h Gn = 


CGnkGX „ D 0.8 


0.3 


dGn 


ReGX n p rGX n 


(5) 


where Cg„ the gas-side correlation coefficient given as input and the Reynolds number is 
defined by 


evaluated via 


or 


„ . _ 4 WG P GSn 

CGXn irdGnUGXn T GXn 

(6) 

TGXn = f( P GSn’iGXn) 

(7) 

P GSn = f( P GS n ,iGSn) 

(8) 

trauoier coefficient is determined the wall heat 

flux can be 

<7n = h Gn (TGAW n - TcWn ) 

(9) 

<?n = „ (*GAW n lGW n ) 

^PGXn 

(10) 


The adiabatic wall temperature and gas-side heat transfer coefficient, calculated from equa- 
tions (4) and (5), or wall heat flux calculated using equations (9) and (10) will be used in the 
conduction module to evaluate a revised wall temperature distribution. It should be noted 
that the formulation given by equations (5-10) yields an approximate value for the wall heat 
flux. To obtain a more accurate value for the wall heat flux a boundary layer model should be 
implemented. The procedure for interfacing a boundary layer module to the present model will 
be described later. Next, attention will be focused on calculating the coolant-side properties 
and heat transfer coefficient. 
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For the first station the coolant stagnation enthalpy, static pressure and static density 
are set equal to the stagnation enthalpy, pressure, and density at the entrance to the cooling 
channel (i.e., icOi =, *C0> R CS\ = p CQ and PCS\ = PCO )• F° r the other stations, the coolant 
stagnation enthalpy is calculated via 


iCOn = iC0 n _! + 


( On 1 + <?n- l)AS n -i. n 

2 W C 


( 11 ) 


where AS n _i, n is the distance between two neighboring stations n — 1 and n and is the 
heat transferred per unit length of the cooling channel from the hot gases to the coolant at 
station n (calculated from the conduction subroutine at iteration j — 1). For the first iteration 
at station n, q £ -1 in equation (11) is not known; therefore the following equation is used to 
evaluate the stagnation enthalpy 


iCOn = *C0n— 1 + 


?n-lASn-l,n 


( 12 ) 


Note that q n - 1 in equations (11) and (12) are the heat transfer per unit length of cooling 
channel at the previous station. 

The coolant velocity is calculated from the following equation: 


VCSn = 


Wc 

PCSn A C n Nn 


(13) 


Note that pcs n > is set equal to pco n l° r the first station, and for the other stations is evaluated 
using the GASP or WASP programs [3, 7] based on the static pressure and enthalpy at the 
previous iteration, i.e., 


PCSn = P( P CSn ’ l CS n ) ( 14 ) 

At the first iteration, however, it is set equal to the static density of the previous station 

( PCSn ~ PCSn-l )• 

Once the coolant velocity is determined, the static enthalpy can be calculated using the 
following equation: 




ics n = iC 0 „ - 


cs „ 
%9cJ 


(15) 


The coolant static and reference Reynolds numbers, respectively, are given by: 


and 


Recs n = 


Wcdc n 
ACn N nPCS n 


(16) 


RecXn - Re CSn 


f PCWn \ / PCSn \ 

\PCSn J \PCW n ) 


(17) 


where pcs n is a function of Pcs n and ics n and is calculated using the GASP program [3], or 
the WASP program [7] if the coolant is water. Note also that dc n is the coolant hydraulic 
diameter at station n. To employ a better value for the Reynolds number, an average Reynolds 
number between the entrance and exit to each station is evaluated, i.e., 


Re C S Avg . = 0.5 (RecSn + RecSn-i) 


(18) 
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Re C X Avg . = 0.5 (Recx n + RecXn-i) (19) 

The Reynolds number in the cooling channel is within the turbulent flow range; hence, the 
Colebrook equation [10] is used to calculate the friction factor. This equation is given by: 


1 nnl f e 2.5226 \ 

■v/7 ° g \^3.7065£> + Re C x Avg V7 ) 


( 20 ) 


This implicit equation has been shown to be very closely approximated by the explicit formula 

ini 


-L = -2.01og 


e 

3.7065D 


5.0452 ( 1 /eX 1 - 1098 5.8506 \ 

Recx Avg . g y2.8257 UJ + R*°cZl,. ) 


( 21 ) 


The correlation given by equation (21) is only valid for straight channels. To include the 
curvature effect, the friction factor obtained from equation (21) must be multiplied by the 
curvature factor given by Ito’s correlation [12]: 


<t>Cur . — 

where rc n is the hydraulic radius of cooling channel Rcur. n is the radius of curvature. The 
curvature factor given by equation (22) is valid when Recx Avg ,{ ) 2 > 6, otherwise, <pcur . = 
1. 

Once the friction factors are determined, the viscous pressure drop between stations n — 1 
and n is calculated using Darcy’s law [13] which is given by: 


\Recx 


r C n 


Avg. 


Rcn 


! 4/^ 


( 22 ) 


(apc..,,,)/ = £ ( 7 C : + + Zt ',‘ ) (Vcs ■ + ifcs - )2As - 1 '" 

and the momentum pressure drop is calculated via 


(23) 


/ 2 \ w c ( 1 1 

(APCS.- 1 , JM - + {NAc)n ) 9c \(pCsA C N) n (p C sA C N) n . l 

An average value of variables between stations n and n — 1 are used to improve the accuracy. 
The static pressure at each station is calculated based on the viscous and momentum pressure 
drops and is given by: 



PCSn = P CS n _i - [(APcS„_ lt J/ + (&PCS n -l,n)M 


(25) 


Once the coolant static pressure is determined, the coolant wall properties which are functions 
of the static coolant pressure Pcs n and wall temperature, i.e., 


C pcw n . Pcw n , kcw n , *cw n = f ( Pcs n - T cw„ ) (26) 

are evaluated using the GASP or WASP programs. It should be noted that the wall temper* 
ature is not constant at a given station; hence, three coolant wall properties which are based 
on the lower, upper and side wall temperatures are determined. The reference and adiabatic 
wall enthalpies at the station are, respectively, calculated from the following equations [8] 
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icx n = °- 5 (*CS„ + icw n ) + 0.194(ico n _ ! CS„ ) 


(27) 


and 


iCAWn = icSn + {Prcx) 1/Z (iC0n ~ *CS U ) ( 28 ) 

The adiabatic wall temperature is a function of the coolant static pressure and the adiabatic 
wall enthalpy and is evaluated using the GASP program [3). Note that the Prandtl number in 
equation (26) is expressed by: 


where 


Prcx = 


kcx 


(29) 


C pcx , HCX, k C X = f{PcS,icx) 

Three correlations may be used to evaluate the heat transfer coefficients in the cooling 
channels. The simplest one is given by the following correlation (see [8. 9]): 

Nu - Cc n Re°c 8 x P r cx ( 3 °) 

Most recently, a new correlation is presented in [14, 15). In this correlation the Nusselt 
number is given by: 

P- = Cc.fleOW' 4 (31) 

Nu r 

where 


jVu r = ip 


- 0.55 


'ip — 1 + i{Tw — Ts) 


and 

_ * ^ _ I (§t)p 

7 " pd T P - p(%) T 

Properties for the above correlation are based on the coolant static temperature Tcs , and 
static pressure Pcs ■ Correlations described by equations (30) and (31) give inaccurate results 
when the coolant is liquid oxygen. A correlation, specifically for oxygen has been proposed 
[16]. This correlation is given by: 


Nu = Cc n RecsP r ° A 



where Peri = 731.4 psia is the critical pressure and 


(32) 


__ iCW ~ i£S 

P Tew - Tcs 
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The properties in the above correlations are calculated using the GASP program [3], or the 
WASP program [7] if the coolant is water. It should also be noted that there are three 
coolant heat transfer coefficients and adiabatic wall temperatures. They are for the top, side, 
and bottom walls of the cooling channel. The variable heat transfer coefficient is due to the 
variable wall temperature in the cooling channel. The coolant reference and adiabatic wall 
enthalpies are also functions of wall temperature and are larger for the surface nodes closer 
to the bottom of the cooling channel. The correlation factors for the heat transfer coefficient, 
Cc n , in equations (30) and (31) are usually equal to 0.023 for most coolants. When the coolant 
is liquid oxygen, however, a factor of 0.0025 is used in equation (32). 

The correlations given by equations (30)-(32) are for fully developed turbulent flow in a 
smooth and straight tube (channel). To include the effect of the entrance region, they are 
multiplied by the following coefficient [17]: 


4>Ent. = 2.88 ( Er= ^ - — -) 


- 0.325 


(33) 


Other entrance effect factors for different types of cooling channel entrances reported in [17] 
are given by: 


<t>Ent.= 1 + ( ^Hc 5ll+1 ) 


for a 90° bend entrance. Taylor [18] suggested the following correction factors: 

<t>Ent. = (Tiv/Tft)! 1 - 59 ^ r«i AS M+i/dcn)] 

for straight tube and 

<pEnt. = (TWr/r t )I 1 - M/ ^ta l A5 <li+1 /d Cn )] j + 5/ 
for a 90° bend entrance. The correction factor for the curvature effect is given by [19]: 


(34) 


(35) 


(36) 


< fiCur . — 


Re cx Au9 . 



21 i 1 / 20 


(37) 


where rc n is the hydraulic radius of cooling channel, Rcur. n is the radius of curvature, the sign 
(+) denotes the concave curvature and the sign (-) denotes the convex one. To incorporate 
the effect of surface roughness on the heat transfer coefficient, a simple empirical correlation 
is suggested by Norris [20]: 


Nu / f 

^ u smooth \ ^smooth 

where n = 0.68Pr 0 215 . For /// smoo th > 4.0 Norris finds that the Nusselt number no longer 
increases with increasing roughness. 

Once the heat transfer coefficients and adiabatic wall temperatures for the hot gas and 
coolant are evaluated, a finite difference model is used to re-evaluate the wall temperature 
distribution. This model has been specifically developed for three-dimensional conduction in 
a rocket thrust chamber and nozzle, as shown in Figure 1. Because of the symmetry of the 
configuration, computations are performed for only one cell (see Figure 3). Since no heat is 
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Figure 3. A Half Cooling Channel Cell. 



(NPHIC) 


Figure 4. Finite Difference Grids Superimposed on Half of a Cooling Channel. 
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transferred to the two sides of the cell, they are assumed to be insulated. A finite difference 
grid is superimposed on the aforementioned cell as shown in Figure 4. In this program the 
number of nodes in the radial direction for different layers and in the circumferential direction 
for the land and channel area must be specified. Thus, the grid size can vary from one layer 
to another. Each node is connected to four neighboring nodes at the same station. It also 
exchanges heat with its counterpoints at two neighboring stations (i.e., stations n + 1 and 
n — 1). The finite difference equation for a node located in the middle of a material is given by 


l _ T i+lJ,n/ R l + T ij-l t n/ R 2 + + T iJ+l,n/ Re 


T • = 
where 


R i = 


1 / R\ + l/i ?2 ■+■ l/Rz + I/Ra + 1 / /?5 + 1/Re 
r A<t> 


1 1 

+ 


Ar (A5”~ 1,n + A5 t n j ni " 1 ) k^ 


R 2 = 


Ar 


1 1 

+ 


(r + f )A<t> (AS"' 1 '" + AS™* 1 ) 1,„, 

r A<f> 


R 3 = 


1 1 
+ 


Ar(A5- 1 ' n + AS^* +1 ) 


/?4 = 


Ar 


1 1 

+ 


(r - #)A 4> (A S^ 1 '" + A S™ +1 ) k l ~/ +hn 


Rs = 


AS?’ 


n,n+l 


h3 


l l 

+ 


\kij tTl k i,j,n + 1 


AS",- 1,n 

Re - — - — 


1 1 

+ 


^ij,n — 


2 i 4 ij )Tl -i \kij, n ki t j,n- 1 J 
(rA0Ar )n+l + (rA<*Ar) n 


( 39 ) 


and 


(rAtf>Ar) n + (rA^Ar ) n _i 

— 2 

and l is the Gauss-Siedel iteration index. Note that the above equation is a three-dimensional 
finite difference equation. The Gauss-Siedel iteration, however, is only performed for the 
nodes on the n-th station and Tij^i and are kept constant during this iteration. 

The value of T^n-i in equation (39) is from the recent march and Tij t7l + 1 from the previous 
march. The conductivity in equation (39) is a function of temperature, i.e., k — k(T). Similar 
equations are derived for other nodes (boundary nodes and nodes at the interface between two 
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different materials). It should be noted that at the boundary nodes, depending on the boundary 




Acf> 


(i* j ”1 »n) 



Figure 5. Resistances between a typical interior node and its adjoining nodes. 


conditions, convective and radiative terms also appear in the nodal balance of energy equation. 
For example, for a node at the inner surface of the nozzle the finite difference equation is given 

by 


TU 1 , j , w /^ 1 + + Tjj,n-l/R§ + Qr 


nl~ 1 


nl~l 


T ■ = 

1 iJ,n 


where 


R i = 


1 / R\ -F 1/^2 + 1/ R3 + 1/ R4 + 1/ R§ + l/ii6 


2rAo 


*2 = 


1 Ar (as”- 1 - + As t 7 +l ) V*Si 

Ar 


1 1 
+ 


(r + y)Ao (A S"j 1,n + AS”j” +1 j \ k i,j,n k i,j-l,n > 

2rAd> 


R 3 = 


1 1 
+ 


Ar (AS”' 1 ’” + AS t ”/ +1 ) 
2 


/?4 = 


* + 


hgrAo (ASy + AS’y) 

k i,j,n+l ) 

AS"- 1 ’" 


i?6= 


ui-1 

1 

PT 


+ 


1 


2j4t,j,n-l \ k ij t n k ij,n - 1 


(40) 
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Ai : j,n ~ ' 
Aij'fi — 1 = 


[( r » + x) A *Tr]„ +1 + [(d + x) 


Ar 


[(rp + jf) A^^] n + [(ro + j) 

2 


Note that equation (40) is used when hot-gas-side heat transfer coefficient is known and wall 
heat flux is evaluated based on the temperature difference, i.e., equation (9). When wall heat 
flux ( q n ) is known, equation (40) becomes 


T Ln = [ T l+Ln/ R 1 + + Tll^JRz + T iij}n+l /R 5 + T l j i n-i/Re+ 

q n A<t>(AS”J l n + A5fj n+1 )/2 + Q r ] /(1/fli + 1 /R 2 + 1/Rz +.I/JI 5 + l/«6) (41) 


where q n is wall heat flux which can be an input of the program or evaluated using equation 
(10). Q r is the radiative heat transfer term which is evaluated based on the Discrete Exchange 
Factor (DEF) method (21, 22, 23. 24] and is given by: 


Qr = 


A<t>(AS™, 1,w + A 


4n 


( m-h 2 

E- 

1=1 


w Sl E Sl DSiS n + ^ ' wgiE g^GiS-n E Sn 


(42) 


1=1 


E Sn and E 9n are surface and gas emissive powers at stations n and are related to their tem- 
peratures via 


E 


Sn 


= ea 


2?rr T 4 
sin (3 n Sn 


Eg, = AK t ,{\ - uo)air r 2 T^ 


DSiSn and DGiS n are total exchange factors between differential surface and gas elements at 
station / to a surface element at station n. The total exchange factor between two elements 
is defined as the fraction of the radiative energy that is emitted from one element and is 
absorbed by the other element via direct radiation and multiple reflections and scatterings 
from surfaces and gas. Procedures for calculating direct and total exchange factors in rocket 
thrust chambers and nozzles are presented in [23] and [24). The radiative heat transfer term, 
given by equation (42), evaluates the radiative energy coming to a surface node from all parts of 
the engine. This is done by numerical integration of the radiative energy incident on the surface 
at station n that is originated (emitted) from station l. The weight factors w s and w g are used 
for numerical integration of surface and gas radiation along axial direction. If the stations 
are equally spaced then the weight factors are the same as those of trapezoidal or Simpson 
methods. In the present application, however, the stations are more concentrated at the throat 
area and are unevenly spaced. The rectangular numerical integration method is suitable when 
stations are not equally spaced and the weight factors are equal to the width of each station, i.e., 
(AS l n ” 1,ri + A5[j n " > ~ 1 )/2. It should be noted that the evaluation of exchange factors DSkS n and 
DGfcSn involves multiple integration (see [23] and [24]) and requires significant computer time. 
Values of these exchange factors depend on the geometry of engine and radiative properties 
of combustion gases. Hence, they can be evaluated using external modules and the resulting 
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exchange factors stored in files for different engines. These files can then be used as inputs to the 
RTE. A separate computer program, namely RTEJDEF (Rocket Thermal Evaluation-Discrete 
Exchange Factor), has been developed for evaluation of the total exchange factors. Note 
that the combustion properties code given by Gordon and McBride [2] does not provide the 
radiative properties of combustion gases. These properties may be obtained from [25] and [26]. 
For example, if the fuel is RP-1, the combustion gas species mole fractions are obtained from 
the combustion code [5], containing 17%C02, 30%CO, 33%H20. 6%OH, 2.5%02, 3%H, 7%H2 
and 1.5%0. Using an integrated average value of the absorption coefficients of these species, 
the overall absorption coefficient is found to be K a = 2.5m” 1 . 

Based on the revised wall temperature, new hot-gas and coolant wall properties, heat 
transfer coefficients and adiabatic wall temperatures are calculated using equations (1) through 
(42). Again, a new wall temperature distribution based on the most recent heat transfer 
coefficients and adiabatic wall temperatures is calculated using the finite difference subroutine 
for heat conduction within the wall. This procedure is repeated until the relative difference 
between the temperature distribution of two consecutive iterations becomes negligibly small. 
After the results for station n converge, the coolant Mach number and entropy as functions 
of static pressure and enthalpy (Mc n ,s Cn = f(Pcs n ^CS n ) ) are evaluated using the GASP 
or WASP programs. Next, the coolant stagnation pressure is evaluated based on the coolant 
entropy and stagnation enthalpy, i.e., Pco n = P(*CO n , scj- The GASP and WASP programs 
do not have explicit expressions for pressure in terms of entropy and enthalpy. Thus, an 
implicit relation for stagnation pressure (i.e., sc n = s (Pco n ,ic n )) with the secant method for 
solving nonlinear equations is used to determine Pco n - the secant method, two initial guesses 
for the stagnation pressures were made (Pi = Pco n _i + 20 and P 2 = Pco n _i - 20) and the 
corresponding entropies si and S2 were determined. The secant method's iterative equation is 
given by: 

n r-> Pfc-l — P/c f AO\ 

Pfc+l = Pk - s k (43) 

s k — 1 s k 

where k is the iteration index. When equation (43) converges (the convergence criterion is 
| Sk — sc n \ < 0.0001), the coolant stagnation is set equal to the latest value of P*.. Finally, the 
coolant stagnation temperature is determined based on the coolant stagnation pressure and 
enthalpy (Tco n = T(Pco n ^co n ))- 

The program then marches axially and performs similar calculations (i.e., equations (1) 
through (43)) for all stations. Once the results of the last station (station m) converged, the 
results of this march are compared to those of the previous march. If the relative differences 
between the results of two consecutive marches is less than the axial convergence criterion the 
program stops, otherwise it continues its axial marches until convergence is achieved. The 
effect of axial conduction can be eliminated by setting the axial convergence criterion greater 
than one or setting the maximum number of passes equal to one. A complete flow chart of 
RTE is presented in [1]. 

The method described here, i.e., axial marches along axial direction, has several advan- 
tages over the direct solution of a three-dimensional finite difference formulation. First, it 
converges very quickly. Second, it requires less memory. Third, it allows the user to control 
the importance of axial conduction by allowing for different convergence criterion between 
the axial and radial and circumferential directions. For example, in analysis of a thin-walled, 
radiatively-cooled, low-pressure engine, axial conduction is negligible. In this case one might 
set the convergence accuracy to 5% in the axial direction and 0.1% in the other directions. In 
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the case of a thick-walled, regeneratively-cooled, high-pressure engine, axial conduction may 
be significant. Thus, the accuracy in the axial direction may be set to 0.1% and 0.1% in the 
other directions. 


Results of a Typical Run 

RTE is used to determine the temperature distribution and heat transfer characteristics of 
a Liquid Oxygen/ Liquid Hydrogen rocket engine. The engine has the following specifications: 


Fuel 

lh 2 

Oxidant 

lo 2 

Coolant 

lh 2 

Chamber stagnation pressure Pgo 

1600 psi 

Coolant stagnation pressure Pco 

2400 psi 

Fuel flow rate 

35.412 lb/sec 

Coolant flow rate 

5.059 lb /sec 

Fuel/Oxidant Mixture ratio 

5.9957 

Coolant stagnation temperature 

110 R 

Number of cooling channels 

100 


The engine is subdivided into 29 stations. Table 1 shows dimensions of the engine and 
some thermal characteristics at each station (see Figure 3 for notation). Note that dimensions 
given in Table 1 are in inches. Also, DCIN = 0.035 in. remains constant along the engine. 
The outer surface radiates to space and its emissivity is 0.9. The thermal conductivities of 
wall materials, i.e., nickel and copper are functions of temperature. 

The resulting wall temperature distributions for stations 1, 9, 16 and 29 are shown in 
Figure 6. A close examination of the temperature distributions reveals that the temperature 
gradient is relatively large in radial direction , especially for station 9. This may also be true 
for any other high pressure thrust chamber. Also, the results shown in Figure 6 can be used 
to optimize the cooling channel aspect ratio. For example, there is no temperature gradient at 
the upper section of the cooling channel in Figure 6a (station 9, throat). Hence, the cooling 
channel can be shorten slightly without changing the overall heat transfer to the coolant. 


HOT-GAS-SIDE BOUNDARY LAYER ANALYSIS INTERFACE 

The convective heat transfer coefficients and heat fluxes for the hot-gas-side of the RTE are 
evaluated based on a tube-like correlation [8], see equation (7). To obtain more accurate results. 
RTE can be linked to a nozzle flow and boundary layer analysis program. The procedure for 
linking RTE to TDK (Two-Dimensional Kinetics Nozzle Performance Computer Program [5]) 
is described in this section. A similar approach may be implemented to link RTE to other 
nozzle boundary layer analysis programs. 

The flowchart for the iterative procedure for linking RTE to TDK is shown in Figure 7. In 
this approach, initially, the wall fluxes and temperatures are evaluated by running RTE under 
an unkn own wall heat flux condition. In this run, RTE uses it internal hot-gas-side routines. 
The wall temperatures calculated by RTE are then used in the inputs of TDK. Using one 
of TDK’s boundary layer modules (BLM or MABL)[5], a new wall heat flux distribution is 
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evaluated. The wall heat flux distribution is inserted into the RTE inputs. This time, since 
the hot-gas-side heat fluxes are known, RTE bypasses all hot-gas-side calculations (e.g., its 
CET subroutine and hot-gas-side heat transfer coefficient correlations) and calculates the wall 
temperature distribution. The new wall temperature distribution along the axial direction is 
then input to TDK and a new heat flux distribution is calculated. This iterative procedure 
continues until convergence is reached. 



Figure 7. Flowchart of RTE-TDK Interface. 


The RTE-TDK model is used to predict wall heat fluxes and temperatures of the LO/LH 
engine presented in the previous section. The resulting wall heat flux and temperature distri- 
butions for both RTE and RTE-TDK calculations are shown in Figures 8 and 9. As shown 
in these figures, the heat flux and temperature distribution when the boundary layer module 
is used are consistently below those calculated via hot-gas-side heat transfer coefficient, i.e., 
equation (5). The reduction of heat flux and temperature is due to the relaminarization of 
accelerating flow. 
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Figure 8. Wall heat fluxes for a H2 cooled engine based on RTE and RTE-TDK models. 



Figure 9. Wall temperatures for a H2 cooled engine based on RTE and RTE-TDK models. 



Table 1: Parameters of thrust chamber and nozzle at different stations. 


Station 

X 

DG 

CCW 

CCH 

THKNS 

1 

3.55 

6.978 

0.05 

0.133 

0.368 

2 

2.75 

5.898 

0.05 

0.123 

0.358 

3 

2.0 

4.886 

0.05 

0.113 

0.348 

4 

1.4 

4.076 

0.05 

0.104 

0.339 

5 

0.9 

3.402 

0.05 

0.1 

0.335 

6 

0.559 

2.942 

0.04 

0.1 

0.335 

7 

0.3 

2.692 

0.04 

0.1 

0.335 

8 

0.1 

2.610 

0.04 

0.1 

0.335 

9 

0 . 

2.6 

0.04 

0.1 

0.335 

10 

-0.1 

2.608 

0.04 

0.1 

0.335 

11 

-0.274 

2.656 

0.04 

0.1 

0.335 

12 

-0.506 

2.746 

0.04 

0.1 

0.335 

13 

-0.906 

3.924 

0.05 

0.1 

0.335 

14 

-1.306 

3.092 

0.05 

0.1 

0.335 

15 

-1.706 

3.264 

0.05 

0.104 

0.339 

16 

-1.906 

3.344 

0.05 

0.113 

0.348 

17 

-2.106 

3.432 

0.05 

0.123 

0.358 

18 

-2.306 

3.516 

0.05 

0.125 

0.36 

19 

-2.506 

3.602 

0.05 

0.125 

0.36 

20 

-2.906 

3.77 

0.05 

0.125 

0.36 

21 

-3.306 

3.94 

0.05 

0.125 

0.36 

22 

-4.106 

4.236 

0.05 

0.125 

0.36 

23 

-4.506 

4.358 

0.05 

0.125 

0.36 

24 

-5.506 

4.6 

0.05 

0.125 

0.36 

25 

-6.506 

4.744 

0.05 

0.125 

0.36 

26 

-7.572 

4.8 

0.05 

0.125 

0.36 

27 

-8.35 

4.8 

0.05 

0.125 

0.36 

28 

-9.0 

4.8 

0.05 

0.125 

0.36 

29 

-9.875 

4.8 

0.05 

0.125 

0.36 


CONCLUDING REMARKS 

The numerical model for a rocket thermal analysis code (RTE) has been discussed. This 
model allows temperature variation along three directions: axial, radial and circumferential. 
The numerical results presented show that there is a large temperature gradient in the axial 
direction for engines with a high chamber pressure. The resulting wall temperature distribution 
can be used to optimize the cooling channel aspect ratios 

The RTE needs to be modified further to incorporate a wide range of cooling channel shapes 
and a CFD model for the cooling channel flow analysis. Efforts are presently under way to 
include these items in the RTE. 
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i 

NOMENCLATURE 


A 

C 

C p 

d 

DG kS n 
DSkS n 
e 
E 

f 

9c 

h 

i 

J 

k 

K t 

m 

N 

P 

Pr 

9 

Qr 

r 

RCur. 

Rn 

Re 

s 

T 

V 

W 

w 

x 


area 

correlation factor for heat transfer coefficient 

specific heat 

diameter 

total exchange factor between gas and surface differential elements 
total exchange factor between two surface differential elements 
cooling channel surface roughness 
surface and gas emissive power 
friction factor 

gravitational constant. 32.2 ft.lb m /lby.s 2 

heat transfer coefficient 

enthalpy 

work/heat proportionality factor 
conductivity 

total extinction coefficient 

total number of axial stations 

total number of cooling channels 

pressure 

Prandtl number 

heat flux 

radiative heat transfer at inner surface 
radius 

radius of curvature 

thermal resistance 

Reynolds number 

entropy 

temperature 

velocity 

weight flow 

weight factor for discrete exchange factor method 
station position in longitudinal direction 


Greek Symbols 

3 angle between a vector normal to the nozzle surface and axial direction 

/\S length of cooling channel between two stations 

Ap pressure drop 

Ar radial mesh size 

Acf) circumferential mesh size 

e convergence criteria or error limit 

^ dynamic viscosity 
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I 


p 

a 

4 > 


density 

Stefan-Boltzmann coefficient 

entrance and curvature effect correction factors 


Subscripts 


A 

Avg. 

C 

Cur. 

f 

G 

i 

j 

k 

M 

n 

r 

S 

s 

W 

X 

0 


adiabatic 

average 

coolant 

curvature 

viscous or friction 

gas 

node i 

node j 

secant method iteration number 

momentum 

related to station n 

radiation 

static 

surface 

wall 

reference 

stagnation 


Superscripts 

j iteration number 

l iteration number for conduction model 

n related to station n 
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